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ABSTRACT 



Context. Current models of the size- and radial evolution of dust in protoplanetary disks generally oversimplify either the radial evo- 
lution of the disk (by focussing at one single radius or by using steady state disk models) or they assume particle growth to proceed 
monodispersely or without fragmentation. Further studies of protoplanetary disks - such as observations, disk chemistry and structure 
calculations or planet population synthesis models - depend on the distribution of dust as a function of grain size and radial position 
in the disk. 

Aims. We attempt to improve upon current models to be able to investigate how the initial conditions, the build-up phase, and the 
evolution of the protoplanetary disk influence growth and transport of dust. 

Methods. We introduce a new version of the model of Brauer et al. (2008) in which we now include the time-dependent viscous 
evolution of the gas disk, and in which more advanced input physics and numerical integration methods are implemented. 
Results. We show that grain properties, the gas pressure gradient, and the amount of turbulence are much more influencing the evolu- 
tion of dust than the initial conditions or the build-up phase of the protoplanetary disk. We quantify which conditions or environments 
are favorable for growth beyond the meter size barrier. High gas surface densities or zonal flows may help to overcome the problem 
of radial drift, however already a small amount of turbulence poses a much stronger obstacle for grain growth. 
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The question of how planets form is one of the key questions 
in modern astronomy today. While it has been studied for cen- 
turies, the problem is still far from being solved. The agglom- 
eration of small dust particles to larger ones is believed to be at 
least the first stage of planet formation. Both laboratory exper- 
iments (Blum et al. 2000) and observations of the 10 //m spec- 
tral region (Bouwman et al. 2001; van Boekel et al. 2003) con- 
clude that grain growth must take place in circumstellar disks. 
The growth from sub-micron size particles to many thousand 
kilometer sized planets covers 13 orders of magnitude in spatial 
scale and over 40 orders of magnitude in mass. To assemble a 
single 1 km diameter planetesimal requires the agglomeration of 
about 10 27 dust particles. These dynamic ranges are so phenom- 
enal that one has to resort to special methods to study the growth 
of particles though aggregation in the context of planet(esimal) 
formation. 

A commonly used method for this purpose makes use of 
particle size distribution functions. The time dependent evo- 
lution of these particle size distribution functions has been 
studied by Weidenschilling (1980), Nakagawa et al. (1981), 
Dullemond & Dominik (2005), Brauer et al. (2008a) (hereafter 
B08a) and others. It was concluded that dust growth by co- 
agulation can be very quick initially (in the order of thousand 
years to grow to centimeter sized aggregates at 1 AU in the so- 
lar nebula), but it stalls around decimeter to meter size due to 
what is known as the "meter size barrier": a size range within 
which particles achieve large enough velocities to undergo de- 
structive collisions and fast radial inward drift toward the central 
star (Weidenschilling 1977; Nakagawa et al. 1986). 
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Wbile the existence of this meter size barrier (ranging in fact 
from a couple of centimeters to tens of meters at 1 AU) has been 
known for over 30 years, a thorough study of this barrier, includ- 
ing all known mechanisms that induce motions of dust grains 
in protoplanetary disks, and at all regions in the disk, for vari- 
ous conditions in the disk and for different properties of the dust 
(such as sticking efficiency and critical fragmentation velocity), 
has been only undertaken recently (see B08a). It was concluded 
that the barrier is indeed a very strong limiting factor in the for- 
mation of planetesimals from dust, and that special particle trap- 
ping mechanisms are likely necessary to overcome the barrier. 

However, this work was based on a static, non-evolving 
gas disk model. It is known that over the duration of the 
planet formation process the disk itself also evolves dra- 
matically (Lynden-Bell & Pringle 1974; Hartmann et al. 1998; 
Hueso & Guillot 2005), which may influence the processes of 
dust coagulation and fragmentation. It is the goal of this paper to 
introduce a combined disk-evolution and dust-evolution model 
which also includes additional physics: we include relative az- 
imuthal velocities, radial dependence of fragmentation critical 
velocities and the Stokes-drag regime for small Reynolds num- 
bers. 

The aim is to find out what the effect of disk formation and 
evolution is on the process of dust growth, how the initial con- 
ditions affect the final outcome, and whether certain observable 
signatures of the disk (for instance, its degree of dustiness at 
a given time) can tell us something about the physics of dust 
growth. 

Moreover, this model will serve as a supporting model for 
complementary modeling efforts such as the modeling of ra- 
diative transfer in protoplanetary disks (which needs informa- 
tion about the dust properties for the opacities) and modeling of 
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the chemistry in disks (which needs information about the total 
amount of dust surface area available for surface chemistry). In 
this paper we describe our model in quite some detail, and thus 
provide a basis for future work that will be based on this model. 

Furthermore, additional physics, such photoevaporation or 
layered accretion can be easily included, which we aim to do in 
the near future. 

As outlined above, this model includes already many pro- 
cesses which influence the evolution of the dust and the gas 
disk. However, there are several aspects we do not include 
such as back-reactions by the dust through opacity or collec- 
tive effects Weidenschilling (1997), porosity effects (Ormel et al. 
2007), grain charging (Okuzumi 2009) or the "bouncing barrier" 
(Zsom et al., in press). 

This paper is organized as follows: Section 2 will describe 
all components of the model including the radial evolution of 
gas and dust, as well as the temperature and vertical structure of 
the disk and the physics of grain growth and fragmentation. In 
Section 3 we will compare the results of our simulations with 
previous steady-state disk simulations and review the aforemen- 
tioned growth barriers. As an application, we show how different 
material properties inside and outside the snow line can cause 
a strong enhancement of dust within the snow line. Section 4 
summarizes our findings. A detailed description of the numeri- 
cal method along with results for selected test cases can be found 
in the Appendix. 

2. Model 

The model presented in this work combines a ID viscous gas 
disk evolution code and a dust evolution code, taking effects of 
radial drift, turbulent mixing, coagulation and fragmentation of 
the dust into account. We model the evolution of gas and dust 
in a vertically integrated way. The gas disk is viscously evolv- 
ing after being built up by in-falling material from a collapsing 
molecular cloud core. 

The radial distribution of grains is subject to gas drag, radial 
drift, and turbulent mixing. To which extend each effect con- 
tributes, depends on the grain/gas coupling of each grain size. By 
simultaneously modeling about 100-200 different grain sizes, 
we are able to follow the detailed evolution of the dust sub-disk 
being the superposition of all sizes of grain distributions. 

So far, the evolution of the dust distribution depends on the 
evolving gas disk but not vice versa. A completely self consis- 
tently coupled code is a conceptually and numerically challeng- 
ing task which will be the subject of future work. 

2.1. Evolution of gas surface density 

Our description of the viscous evolution of the gas disk follows 
closely the models described by Nakamoto & Nakagawa (1994) 
and Hueso & Guillot (2005). In this paper we shall therefore be 
relatively brief and put emphasis on differences between those 
models and ours. 

The viscous evolution of the gas disk can be described by the 
continuity equation, 
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where the gas radial velocity u g is given by (see 
Lynden-Bell & Pringle 1974) 
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u„ = - 
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Xco 
oo p g (z)dz is the gas surface density, r the radius along 

the disk mid-plane and v„ the gas viscosity. The source term on 
the right hand side of Eq. 1, denoted by S g can be either infall of 
material onto the disk or outflow. 

The molecular viscosity of the gas is too small to account for 
relevant accretion onto the star, the time scale of viscous evolu- 
tion would be in the order of several billion years. Observed ac- 
cretion rates and disk lifetimes can only be explained if turbulent 
viscosity drives the evolution of circumstellar disks. Therefore 
Shakura & Sunyaev (1973) parameterized the unknown viscos- 
ity as the product of a velocity scale and a length scale. The 
largest reasonable values for these scales in the disk are the pres- 
sure scale height H v 

c s 



(3) 



and the sound speed c s . Therefore the viscosity is written as 

y g = ffCs Hp, 



(4) 



where a is the turbulence parameter and a < 1 . 

Today it is generally believed that disks transport angu- 
lar momentum by turbulence, however the source of this tur- 
bulence is still debated. The magneto-rotational instability is 
the most commonly accepted candidate for source of turbu- 
lence (Balbus & Hawley 1991). Values of a are expected to 
be in the range of 10~ 3 to some 10 -2 (see Johansen & Klahr 
2005; Dzyurkevich et al., in press). Observations confirm 
this range with higher probability for larger values (see 
Andrews & Williams 2007). 

If the disk becomes gravitationally unstable, large scale 
mechanisms of angular momentum transport such as through 
the formation of spiral arms come into play. The stability of the 
disk can be described in terms of the Toomre parameter (Toomre 
1964) 

G = -^-- (5) 
nG £g 

Values below a critical value of Q cr = 2 describe a weakly unsta- 
ble disk, which forms non-axisymmetric instabilities like spiral 
arms. Q values below unity lead to fragmentation of the disk. 

The effect of these non-axisymmetric structures is to trans- 
port angular momentum outward and rearranging the surface 
density in the disk so as to counteract the unstable configuration. 
This mechanism is therefore to a certain extent self-limiting. 
Values above Q a are not influenced by instabilities, values be- 
low Q a form instabilities which increase Q until the disk is 
marginally stable again. Our model approximates this mecha- 
nism by increasing the turbulence parameter a according to the 
recipe of Armitage et al. (2001), 



a{r) = a + 0.01 



Qc 



min«2(r), Q CI ) 



(6) 



(2) 



where ao is a free parameter of the model which is taken to be 
10 -3 , unless otherwise noted. 

During the time of infall onto the disk, we use a constant, 
high value of cr = 0.1 to mimic the effective redistribution of 
surface density during the infall phase which also increases the 
overall stability of the disk. Once the infall stops, we gradually 
decrease the turbulence parameter on a timescale of 10 000 years 
until it reaches its input value. 

Eq. 1 is a diffusion equation, which is stiff. This means, one 
faces the problem that the numerical step of an explicit integra- 
tion scheme goes oc Ar 2 (where Ar is the radial grid step size) 
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which would make the simulation very slow. One possible solu- 
tion to this problem is using the method of implicit integration. 
This scheme keeps the small time scales of diffusion i.e. the fast 
modes in check. We are not interested in these high frequency 
modes, but they would become unstable if we used a large time 
step. With an implicit integration scheme (see Section A.l) the 
time step can be chosen larger without causing numerical insta- 
bilities, thus increasing the speed of the computation. 



2.2. Radial evolution of dust 

If the average dust-to-gas ratio in protoplanetary disks is in the 
order of 1CT 2 (as found in the ISM), then the dust does not dy- 
namically influence the gas, while the gas strongly affects the 
dynamics of the dust. 

Thermally, however, the dust has potentially a massive influ- 
ence on the gas disk evolution through its opacity, but we will not 
include this in this paper. Therefore the evolution of the gas disk 
can, in our approximation, be done without knowledge of the 
dust evolution, while the dust evolution itself does need knowl- 
edge of the gas evolution. 

There might be regions, where dust accumulates (such as the 
mid-plane of the disk or dead-zones and snow-lines) and its in- 
fluence becomes significant or even dominant but we will not 
include feedback of such dust enhancements onto the disk evo- 
lution in this paper. 

For now, we want to focus on the equations of motion of dust 
particles under the assumption that gas is the dominant material 
by mass. The interplay between dust and gas can then be de- 
scribed in terms of a dimensionless coupling constant, the Stokes 
number which is defined as 



St 



Ted' 



(7) 



where T e d is the eddy turn-overtime and r s is the stopping time. 

The stopping time of a particle is defined as the ratio of 
the momentum of a particle divided by the drag force act- 
ing on it. There are four different regimes for the drag force 
which determine the dust-to-gas coupling (see Whipple 1972; 
Weidenschilling 1977). Which regime applies to a certain parti- 
cle, depends on the ratio between mean free path A m f p of the gas 
molecules and the dust particle size a (i.e. the Knudsen number) 
and also on the particle Reynolds-number Re = 2au/v mo \ with 
v mo i being the gas molecular viscosity 



v mol - - U A m f p , 



(8) 



u the mean thermal velocity. The mean free path is taken to be 



tmfp 



(9) 



where n denotes the mid-plane number density and <x H2 = 2 x 
1(T 15 cm 2 . 



The different regimes 1 are 



T s = 



Pg« 
9v m ,,i p. 



for A m f p /a > | Epstein regime 



for Re < 1 



Stokes regime 1 



Jo.6 ^1° „o.4 for 1 < Re < 800 Stokes regime 2 



(10) 



^ for Re > 800 Stokes regime 3 

p g « 



Here u denotes the velocity of the dust particle with respect to 
the gas, u = c s yfnJE denotes the mean thermal velocity of the 
gas molecules, p s is the solid density of the particles and p g is 
the local gas density. 

For now, we will focus on the Epstein regime. To calculate 
the Stokes number, we need to know the eddy turn-over time. 
As noted before, our description of viscosity comes from a di- 
mensional analysis. We use a characteristic length scale L c and a 
characteristic velocity scale V c of the eddies. This prescription is 
ambiguous in a sense that it does not specify if angular momen- 
tum is transported by large, slow moving eddies or by small, fast 
moving eddies, that is 



(a"V c ) ■ (ar^Lc). 



(ID 



This is rather irrelevant for the viscous evolution of the gas 
disk, since all values of q lead to the same viscosity, but if we 
calculate the turn-over-eddy time, we get 



Ted 



2nL c 



1 



(12) 



The Stokes number and therefore the dust-to-gas coupling as 
well as the relative particle velocities strongly depend on the 
eddy turnover time and therefore on q . In this work q is taken to 
be 0.5 (following Cuzzi et al. 2001; Schriipler & Henning 2004) 
which leads to a turn-over-eddy time of 



The Stokes number then becomes 
St = ^— ■ - 



for a < -Arafp. 



(13) 



(14) 



The overall radial movement of dust surface density Ed can now 
be described by an advection-diffusion equation, 



(15) 



where the total Flux, F M has contributions from a diffusive and 
an advective flux. The diffusive part comes from the fact that 
the gas is turbulent and the dust couples to the gas. The dust 
is therefore turbulently mixed by the gas. Mixing counteracts 
gradients in concentration, in this case it is the dust-to-gas ratio 
of each size that is being smoothed out by the turbulence. The 
diffusive flux can therefore be written as 



iff 



n 9 i U 



(16) 



1 It should be noted that "Stokes regime" refers to the regime where 
the drag force on a particle is described by the Stokes law - this is not 
directly related to the Stokes number. 
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The ratio of gas diffusivity D g over dust diffusivity is called 
the Schmidt number. We follow Youdin & Lithwick (2007), who 
derived 

Sc= — = 1 + St 2 , (17) 
Dd 

and assume the gas diffusivity to be equal to the turbulent gas 
viscosity v g . 

The second contribution to the dust flux is the advective flux, 



adv 



£ d • Mr 



(18) 



where u r is the radial velocity of the dust. There are two contri 
butions to it, 

Kg 2m„ 



i + st 2 st + sr 



(19) 



The first term is a drag term which comes from the radial move- 
ment of the gas which moves with a radial velocity of « a , given 
by Eq. 2. Since the dust is coupled to the gas to a certain extend, 
the radially moving gas is able to partially drag the dust along. 

The second term in Eq. 19 is the radial drift velocity with 
respect to the gas. The gas in a Keplerian disk does in fact move 
sub-keplerian, since it feels the force of its own pressure gradient 
which is usually pointing inwards. Larger dust grains do not feel 
this pressure and move on a keplerian orbit. Therefore, from a 
point of view of a (larger) dust particle, there exists a constant 
head wind, which causes the particle to loose angular momentum 
and to drift inwards. This depends on the coupling between the 
gas and the particle and is described by the second term in Eq. 
19. u„ denotes the maximum drift velocity of a particle, 



-E A 



dr 



2p g Q. k 



(20) 



which has been derived by Weidenschilling (1977). Here, we 
introduced a radial drift efficiency parameter £<j. This parame- 
ter describes how efficient the radial drift actually is, as several 
mechanisms such as zonal or meridional flows might slow down 
radial drift. This will be investigated in Section 3.5. 

Putting all together, the time dependent equation for the sur- 
face density of one dust species of mass m-, is given by 



<9E 



\_d_ 
dt r dr 



3] 



(21) 



where we have included a source term S d on the right hand side 
which can be positive in the case of infall or re-condensation of 
grains and negative in the case of dust evaporation or outflows. 
This source term does not include the sources caused by coagula- 
tion and fragmentation processes (see Section 2.6. 1). All sources 
will be combined into one equation later which is implicitly in- 
tegrated in an un-split scheme (see Section A. 3). 

Note that both, the diffusion coefficient and the radial veloc- 
ity depend on the Stokes number and therefore on the size of the 
particle. 
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Fig. 1. Total amount of in-fallen surface density as function of 
radius according to the Shu-Ulrich infall model (see Section 2.5) 
assuming a centrifugal radius of 8 AU (solid line), 33 AU 
(dashed line), and 100 AU (dash-dotted line). 



shorter than the radial drift time scale of the gas. The isother- 
mal vertical density structure is then given by 



Pg(z) = Po exp 



1 z 2 ^ 



2 Hi) 



where 



Po = 



(22) 



(23) 



V2^// p 

The viscous heating is given by Nakamoto & Nakagawa (1994) 



g+ =£ g v g r 



~~dr~ 



(24) 



were v g denotes the turbulent gas viscosity and the Kepler 
frequency. 

Nakamoto & Nakagawa (1994) give a solution for the mid- 
plane temperature with an optically thick and an optical thin con- 
tribution, 



(25) 



where we used v g = ac s H p and the approximation tr/p = 
KR/p£g/2. kr and Kp are Rosseland and Planck mean opacities 
which will be discussed in the next section. 

Tin- contains contributions due to stellar or external irradia- 
tion. Here, we use a formula derived by Ruden & Pollack (see 
Ruden & Pollack 1991, App. B), 



3tt I r ) + 2 \ r ) \ r 



d\nH p 
dlnr 



-1 



, (26) 



2.3. Temperature and vertical structure 

The vertical structure can be assumed as being in hydro- 
static equilibrium at all times if the disk is geometrically thin 
(H p (r)/r <k r) and the vertical sound crossing time is much 



with a fixed dln// p /dlnr = 9/7, following Hueso & Guillot 
(2005). 

The main source of opacity is the dust. Due to viscous heat- 
ing, the temperature will increase with surface density. If the 
temperature rises above ~ 1500 K, the dust (i.e. the source of 



T. Birnstiel et al.: Gas- and dust evolution in protoplanetary disks 



5 



opacity) will evaporate, which stops the disk from further heat- 
ing until all dust is vaporized. To simulate this behavior in our 
model, we calculate a gas temperature (assuming a small, con- 
stant value for gas opacity) in the case where the dust temper- 
ature rises above the evaporation temperature. Then r m jd is ap- 
proximated by 



(27) 



only if r m ;d from Eq. 25 would be larger than r evap . 

This is a thermostat effect: if T rises above 1500 K, dust will 
evaporate, the opacity will drop and the temperature stabilizes at 
T - 1500 K. Once even the very small gas opacity is enough 
to get temperatures > 1500 K, all the dust is evaporated and the 
temperature rises further. 

2.4. Opacity 

In the calculation of the mid-plane temperature we use 
Rosseland and Planck mean opacities which are being calculated 
from a given frequency dependent opacity table. The results are 
stored in a look up table and interpolated during the calculations. 
The opacity table is for a mixture of 50% silicates and 50% car- 
bonaceous grains. 

Since these are dust opacities, we convert them from cm 2 /g 
dust to cm 2 /g gas by multiplying the values with the dust-to-gas 
ratio eo, which is a fixed parameter in our model, taken to be 
the canonical value of 0.01. This assumes that the mean opac- 
ity of the gas is very small and that the dust-to-gas ratio does 
not change with time. To calculate the opacity self-consistently, 
the total mass of dust and the distribution of grain sizes has to be 
taken into account, meaning that the dust evolution has a back re- 
action on the gas by determining the opacity. For now, our model 
does not take back-reactions from dust to gas evolution into ac- 
count. Only in the case where the temperature rises above 1500 
K, the drop of opacity due to dust evaporation is considered, as 
described above. 



2.5. Initial infall phase: cloud collapse 

The evolution of the protoplanetary disk also depends on the ini- 
tial infall phase which builds up the disk from the collapse of 
a cold molecular cloud core. This process is still not well un- 
derstood. First similarity solutions for a collapsing sphere were 
calculated by Larson (1969) and Penston (1969). Shu (1977) 
found a self similar solution for a singular isothermal sphere. 
The Larson & Penston solution predicts much larger infall rates 
compared to the inside-out collapse of Shu (mj n m 47 c 3 /G and 
0.975 c^/G respectively). 

More recent work has shown that the infall rates are not con- 
stant over time, but develop a peak of high accretion rates and 
drop to smaller accretion rates at later times. The maximum ac- 
cretion rate is about 13 c^/G if opacity effects are included (see 
Larson 2003). Analytical, pressure-free collapse calculations of 
Myers (2005) show similar behavior but with a smaller peak ac- 
cretion rate of m m = 7.07c 3 /G. They also argue that the effects of 
pressure and magnetic fields will further increase the time scales 
of cloud collapse. 

This initial infall phase is important since it provides the ini- 
tial condition of the disk and also influences the whole simu- 
lation by providing a source of small grains and gas at larger 
distances to the star during later times of evolution. 

It should be noted that several groups perform 3D hydro- 
dynamic simulations of collapsing cloud cores which show 
more complicated evolution (e.g., Banerjee & Pudritz 2006; 



Whitehouse & Bate 2006). However, to be able to study gen- 
eral trends of the infall phase, we use the Shu-model since it is 
adjustable by a few parameters whose influences are easy to un- 
derstand. In this model the collapse proceeds with an infall rate 
of m m = 0.975 c 3 /G which stays constant throughout the col- 
lapse. 

We assume the singular isothermal sphere of mass M c i oul j to 
be in solid body rotation Q. s . If in-falling material is conserving 
its angular momentum, all matter from a shell of radius r s will 
fall onto the star and disk system (with mass m cent (f)) within the 
so called centrifugal radius, 



TcentrW 



G »l cent (f) ' 



(28) 



where G is the gravitational constant and r s = 0.975 ■ c s t/2. The 
path of every parcel of gas can then be described by a ballistic 
orbit until it reaches the equatorial plane. The resulting flow onto 
the disk is 

i i (r,t) = 2pi(r,t)-u z (r,t), (29) 



where 



and 



Gm ceDt (f) 
Uz(r, t) = \\ fi, 



(30) 



p\(r,t) = 



8tt VG m cent (f) r 3 r centr (r, t) fd 2 
as described in Ulrich (1976). Here, fi is given by 



\ (3D 



A* = V 1 - r/r centr (r,t). 



(32) 



The centrifugal radius can therefore be approximated by (cf. 
Hueso & Guillot (2005)) 



r C entri(0 - L4 



10 



■14 , 



3 

( m cent (t) \ I 



3x 10 4 cms- 1 . 



AU. 

(33) 

We admit that this recipe for the formation of a pro- 
toplanetary disk is perhaps oversimplified. Firstly, as shown 
by Visser et al. (2009), the geometrical thickness of the disk 
changes the radial distribution of in-falling matter onto the disk 
surface, because the finite thickness may "capture" an in-falling 
gas parcel before it could reach the midplane. Secondly, star for- 
mation is likely to be messier than our simple single-star axisym- 
metric infall model. And even in such a simplified scenario, the 
Shu inside-out collapse model is often criticized as being unre- 
alistic. However, it would be far beyond the scope of this paper 
to include a better infall model. Here we just want to get a feel- 
ing for the effect of initial conditions on the dust growth, and we 
leave more detailed modeling to future work. 

2.6. Grain growth and fragmentation 

The goal of the model described in this paper is to trace the evo- 
lution of gas and dust during the whole lifetime of a protoplan- 
etary disk, including the grain growth, radial drift and turbulent 
mixing. 

Here, the problem arises that radial drift and coagulation 
"counteract" each other in the regime of St = 1 particles: co- 
agulation of smaller sizes restores the population around St = 1, 
whereas radial drift preferentially removes these particles. To be 
able to properly model this behavior, the time step has to be cho- 
sen small enough if the method of operator splitting is used. 
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Fig. 2. Sources of relative particle velocities considered in this model (Brownian motion velocities are not plotted) at a distance of 
10 AU from the star. The turbulence parameter a in this simulation was 1CT 3 . It should be noted that relative azimuthal velocities 
do not vanish for very large and very small particles. 



The upper limit for this time step can be very small. If larger 
steps are used the solution will "flip-flop" back and forth be- 
tween the two splitted sub-steps of motion and coagulation, and 
the results become unreliable. A method to allow the choice of 
large time steps while preserving the accuracy is to use a non- 
splitted scheme in which the integration is done "implicitly". 
B08a already use this technique to avoid flip-flopping between 
coagulation and fragmentation of grains, and we refer to that pa- 
per for a description of the general method. What is new in the 
current paper is that this implicit integration scheme is extended 
to also include the radial motion of the particles. So now radial 
motion, coagulation and fragmentation are done all within a sin- 
gle implicit integration scheme. See Appendix A. 3 for details. 

2.6.1. Smoluchowski equation 

The dust grain distribution n(m, r, z), which is a function of mass 
m, distance to the star r and height above the mid-plane z, de- 
scribes the number of particles per cm 3 per gram interval in par- 
ticle mass. This means that the total dust density in g cirT 3 is 
given by 

n(m, r,z) ■ m dm. (34) 

o 



With this definition of n(m, r, z), the coagula- 
tion/fragmentation at one point in the disk can be described by a 
general two-body process, 



—n(m,r,z)= II M(m,m',m")x 
at JJo 



(35) 



X n(m , r, z) ■ n(m", r, z) dm' dm", 



where M(m, m', m") is called the kernel. In the case of coagula- 
tion and fragmentation, this is given by 



M(m, m', m") = —K(m', m") ■ 5(m' + m" — m) 

— K(m' , m") ■ 6(m" — m) 

+ —L(m',m") ■ S(m,m',m") 

— L{m' ', m") • 6{m — m"). 



(36) 



For better readability, the dependency of M on radius and height 
above the mid-plane was omitted here. The combined coagula- 
tion/fragmentation kernel consists of four terms containing the 
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coagulation kernel K, the fragmentation kernel L and the distri- 
bution of fragments after a collision S . 

The first two terms in Eq. 36 correspond to gain (masses m! 
and m - m' coagulate) and loss (m coagulates with m') due to 
grain growth. 

The third term describes the fragmentation of masses m and 
m', governed by the fragmentation kernel L and the fourth term 
describes the fact that when masses m' and m" fragment, they 
distribute some of their mass via fragments to smaller sizes. 

The coagulation and fragmentation kernels will be described 
in section 2.6.3, the distribution of fragments, 5, will be de- 
scribed in the next section. 

To be able to trace the size and radial evolution of dust in 
a combined way, we need to express all contributing processes 
in terms of the same quantity. Hence, we will formulate the co- 
agulation/fragmentation equation in a vertically integrated way. 
The vertical integration can be done numerically (as in B08a), 
however coagulation processes are most important near the mid- 
plane, which allows to approximate the kernels as being verti- 
cally constant (using the values at the mid-plane). If the ver- 
tical dust density distribution of each grain size is taken to be 
a Gaussian (see Section 2.6.4, Eq. 51), then the vertical inte- 
gration can be done analytically, as discussed in Appendix A. 2. 
Unlike the steady-state disk models of B08a which have fixed 
surface density and temperature profiles, we need to recompute 
the coagulation and fragmentation kernels (which are functions 
of surface density and temperature) at every time step. Therefore 
this analytical integration also saves significant amounts of com- 
putational time. 

We therefore define the vertically integrated dust surface 
density distribution per logarithmic bin of grain radius, <x(r, a) 



<r(r, a) 



n 

oo 



(r, z, fl) • m • a dz, 



(37) 



where n(r, z, a) and n(r, z, m) are related through m = 4n/3p s a 3 . 
The total dust surface density at any radius is then given by 



Mr) 



-f 

Jo 



<x(r, a) dlno. 



(38) 



Defining cr(r,a) as in Eq. 37 makes it a grid-independent 
density unlike the mass integrated over each numerical bin. This 
way, all plots of cr(r, a) are meaningful without knowledge of 
the size grid that was used. Numerically, however we use the 
discretized values as defined in the appendix. 

In our description of growth and fragmentation of grains, 
we always assume the dust particles to have a constant volume 
density meaning that we do not trace the evolution of porosity 
of the particles as this is currently computationally too expen- 
sive with a statistical code as used in this work. This can be 
achieved with Monte-Carlo methods as in Ormel et al. (2007) or 
Zsom & Dullemond (2008), however these models have do not 
yet include the radial motion of dust and therefore cannot trace 
the global evolution of the dust disk. 

2.6.2. Distribution of fragments 

The distribution of fragments after a collision, S(m,m',m"), is 
commonly described by a power law, 



n(m)dm oc m ^Am. 



(39) 



The value £ has been investigated both experimentally and 
theoretically. Typical values have been found in the range be- 
tween 1 and 2, by both experimental (e.g., Blum & Muench 



1993; Davis & Ryan 1990) and theoretical studies (Ormel et al. 
2009). Unless otherwise noted, we will follow B08a by using the 
value of ij= 1.83. 

In the case where masses of the colliding particles differ by 
orders of magnitude, a complete fragmentation of both parti- 
cles is an unrealistic scenario. More likely, cratering will oc- 
cur (Sirono 2004), meaning that the smaller body will excavate 
a certain amount of mass from the larger one. The amount of 
mass removed from the larger one is parameterized in units of 
the smaller body m s , 

maut=X m $- (40) 

The mass of the smaller particle plus the mass excavated from 
the larger one is then distributed to masses smaller than m s ac- 
cording to the distribution described by Eq. 39. In this work, we 
follow B08a by assuming ;f to be unity, unless otherwise noted. 

2.6.3. Coagulation and fragmentation kernels 

The coagulation kernel K(m\,ni2) can be factorized into three 
parts, 

K(m l ,m 2 ) = Au(m l ,m 2 )o- geo (mi,m2) p c , (41) 



and, similarly, the fragmentation kernel can be written as 
L{m\,mi) — Au{m\,ni2) cr wo (m\,m2) pf. 



(42) 



Here, Au{m\,m2) denotes the relative velocity of the two parti- 
cles, <x ge0 (mi, m2) is the geometrical cross section of the colli- 
sion and p c and pf are the coagulation and fragmentation prob- 
abilities which sum up to unity. In general, all these factors can 
also depend on other material properties such as porosity, how- 
ever we always assume the dust grains to have a volume density 
of p s - 1.6 g cirT 3 . 

The fragmentation probability is still not well known and 
subject of both theoretical (Paszun & Dominik2009; Wada et al. 
2008) and experimental research (see Blum & Wurm 2008; 
Guttler et al. 2009). In this work, we adopt the simple recipe 



Pf = 







1 - 



if Au < Uf — 6u 
if Am > Uf 



(43) 



6u 



else 



with a transition width 6u and the fragmentation speed Uf as free 
parameter which is assumed to be 1 m s -1 , following experi- 
mental work of Blum & Muench (1993) and theoretical studies 
of Leinhardt & Stewart (2009). 

2.6.4. Relative particle velocities 

The different sources of relative velocities considered here are 
Brownian motion, relative radial and azimuthal velocities, tur- 
bulent relative velocities and differential settling. These contri- 
butions will be described in the following, an example of the 
most important contributions is shown in Figure 2. 

Brownian motion, the thermal movement of particles, depen- 
dents on the mass of the particle. Hence, particles of different 
mass have an average velocity relative to each other which is 
given by 



An 



BM 



8k B T{m\ + mj) 
nm\ mi 



(44) 



Particle growth due to Brownian relative motion is most effective 
for small particles. 
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Radial drift, as described in section 2.2 also induces relative 
velocities since particles of different size are differently coupled 
to the gas. The relative velocity is then 



Amrd - \u x (tni) — u r (rti2)\ , 



(45) 



where the radial velocity of the dust, u r is given by Eq. 19. 

Azimuthal relative velocities are induced by gas drag in 
a similar way as radial drift. However while only particles 
(plus/minus 2 orders of magnitude) around St=l are significantly 
drifting, relative azimuthal velocities do not vanish for encoun- 
ters between very large and vary small particles (see Figure 2). 
Consequently, large particles are constantly suffering high ve- 
locity impacts of smaller ones. According to Weidenschilling 
(1977) and Nakagawa et al. (1986), the relative azimuthal veloc- 
ities for gas-dominated drag are 



1 



1 



U + Stf 1 + St 



2 ' 



(46) 



where u n is defined by Eq. 20. 

Turbulent motion as source of relative velocities is discussed 
in detail in Ormel & Cuzzi (2007). They also derive closed form 
expressions for the turbulent relative velocities which we use in 
this work. 

Differential settling is the fifth process we consider con- 
tributing to relative particle velocities. Dullemond & Dominik 
(2004) constructed detailed models of vertical disk structure de- 
scribing the depletion of grains in the upper layers of the disk. 
They show that the equilibrium settling speed for particles in the 
Epstein regime is given by M set t = -z £\ St which can be derived 
by equating the frictional force Fm c = —mu/tfnc an d the vertical 
component of the gravity force, F grav — -mzQj;. To limit the set- 
tling speed to velocities smaller than half the vertically projected 
Kepler velocity, we use 



St, ^ 



(47) 



for calculating the relative velocities. 

Since we do not resolve the detailed vertical distribution of 
particles, we take the scale height of each dust size as average 
height above the mid-plane, which gives 

Akds = \h ■ min(St ; , 1/2) - hj ■ min(St/, 1/2)| ■ Q k . (48) 

The dust scale height is calculated by equating the time scale 
for settling, 

fsett = — (49) 

Msett 

and the time scale for stirring (Dubrulle et al. 1995; 
Schrapler & Henning 2004; Dullemond & Dominik 2004), 



fstir — 



z 

AT 



(50) 



By limiting the vertical settling velocity as in Eq. 47 and by 
constraining the dust scale height to be at most equal to the gas 
scale height, one can derive the dust scale height to be 



H v ■ min 



min(St, 1/2)(1 + St 2 ) 



(51) 



This prescription is only accurate for the dust close to the 
mid-plane, however most of the dust (and hence most of the 
coagulation/fragmentation processes) take place near the mid- 
plane, therefore this approximation is accurate enough for our 
purposes. 



10" 10 ' 
-1.8x10 4 yrs 
-3.0 x 10 4 yrs 
- 1 .0 x 1 6 yrs 




Fig. 3. Evolution of disk surface density distribution (top) and 
midplane temperature (bottom) of the fiducial model described 
in 3.1. 



3. Results 

3. 1 . Viscous evolution of the gas disk 

We will now focus on the evolution of a disk around a T Tauri 
like star. We assume the rotation rate of the parent cloud core 
to be 7 x 10~ 14 s , which corresponds to 0.06 times the break 
up rotation rate of the core. The disk is being built-up from in- 
side out due to the Shu-Ulrich infall model, with the centrifugal 
radius being 8 AU. The parameters of our fiducial model are 
summarized in Table 1 . 

Figure 3 shows how the gas surface density and the mid- 
plane temperature of this model evolve as the disk gets built up, 
viscously spreads and accretes onto the star. It can be seen that 
viscous heating leads to a strong increase of temperature at small 
radii. This effect becomes stronger as the disk surface density 
increases during the infall phase. After the infall has ceased, the 
surface density and therefore also the amount of viscous heating 
falls off. 

This effect also influences the position of the dust evapora- 
tion radius, which is assumed to be at the radius where the dust 
temperature exceeds 1500 K. This position moves outwards dur- 
ing the infall (because of the stronger viscous heating described 
above). Once the infall stops, the evaporation radius moves back 
to smaller radii as the large surface densities are being accreted 
onto the star. 

Figure 4 shows the evolution of accretion rate onto the star, 
stellar mass and disk mass. The infall phase lasts until about 
150000 years. At this point, the disk looses its source of gas 
and small-grained dust and the disk mass drops off quickly until 
the disk has adjusted itself to the new condition. This also ex- 
plains the sharp drop of the accretion rate. The slight increase in 
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0.1 

time [Myr] 



Fig. 4. Evolution of disk mass and stellar mass (solid and dashed 
line on left scale respectively) and accretion rate onto the star 
(dash-dotted line on right scale). Adapted from Figure 5 in 
Hueso & Guillot (2005). 

Table 1. Parameters of the fiducial model. 



parameter 


symbol value 


unit 


turbulence parameter 


a 


io- 3 




irradiation angle 


<P 


0.05 




cloud mass 




0.5 


M Q 


effect, speed of sound in core 


c s 


3 x 10 4 


cm s" 1 


rotation rate of cloud core 


-Q, 


7 x 10~ 14 


s-» 


solid density of dust grains 


Ps 


1.6 


g/cm 3 


stellar radius 


R* 


2.5 


*o 


stellar temperature 




4000 


K 



the accretion rate afterwards comes from the change in a after 
the infall stops (see Section 2.1). Hueso & Guillot (2005) find a 
steeper, power-law decline of the accretion rate after the end of 
the infall phase because their model does not take the effects of 
gravitational instabilities into account. 

3.2. Fiducial model without fragmentation 

We will now focus on the dust evolution of the disk. This fiducial 
simulation includes only grain growth without fragmentation or 
erosion. All other parameters as well as the evolution of the gas 
surface density and mid-plane temperature are the same as dis- 
cussed in the previous section. The evolution of this model is 
visualized in Figure 5. 

The contour levels in Figure 5 show the vertically integrated 
dust surface density distribution per logarithmic bin of grain ra- 
dius, cr(r,a), as defined in Eq. 37. The left column of Figure 5 
shows the results of dust evolution for a steady state (i.e. not 
viscously evolving) gas disk as described in B08a. 

The right column shows the evolution of the dust density dis- 
tribution of the fiducial model, in which the gas disk is gradually 
built up through infall from the parent molecular cloud core, and 
the gas disk viscously spreads and accretes matter onto the star. 
The solid line marks the grain size corresponding to St= 1 at the 
given radius. This grain size will be called ast=i hereafter. In the 
Epstein regime, «st=i is proportional to the gas surface density 
of the disk, which can be seen from Eq. 14. 



There are several differences to the simulations of grain 
growth in steady-state disks, presented in B08a: viscously evolv- 
ing disks accrete onto the star by transporting mass inwards and 
angular momentum outwards. This leads to the fact that some 
gas has to be moving to larger radii because it is 'absorbing' 
the angular momentum of the accreting gas. This can be seen in 
numerical simulations, but also in the self similar solutions of 
Lynden-Bell & Pringle (1974). They show that there is a radius 
R ± which divides the disk between inward and outward moving 
gas. This radius itself is constantly moving outwards, depending 
on the radial profile of the viscosity. 

The radius R+ in the fiducial model was found to move from 
around 20 AU at the end of the infall to 42 AU at 1 Myr and 
is denoted by the dotted line in Figure 5. Important here is that 
small particles are well enough coupled to the gas to be trans- 
ported along with the outward moving gas while larger particles 
decouple from the gas and drift inwards. Those parts of the dust 
distribution which lie below the dotted line in Figure 5 have pos- 
itive flux in the radial direction due to the gas-coupling. 

One might expect that the dotted and dashed lines always co- 
incide for small grains as they are well coupled to the gas, how- 
ever, it can be seen that this is not the case. The reason for this 
is that turbulent mixing also contributes to the total flux of dust 
of each grain size. The smallest particles in between the dotted 
line and the dashed line in the lower two panels of Figure 5 do 
have positive velocities, but due to a gradient in concentration of 
these dust particles, the flux is still negative. 

During the early phases of its evolution, a disk which is built 
up from inside out quickly grows large particles at small radii 
which are lost to the star by radial drift. During further evolution, 
growth timescales become larger and larger (since the dust-to- 
gas ratio is constantly decreasing) while only the small particles 
are mixed out to large radii. 

The radial dependence of the dust-to-gas ratio after 1 Myr is 
shown in Figure 6. These simulations show that the initial con- 
ditions of the stationary disk models (such as shown in B08a 
and in the left column of Figure 5) are too optimistic since they 
assume a constant dust-to-gas ratio at the start of their simula- 
tion throughout the disk which is not possible if the disk is being 
built-up from inside out unless the centrifugal radius is very large 
(in which case the disk would probably fragment gravitationally) 
and grain fragmentation is included to prevent the grains from 
becoming large and start drifting strongly. Removal of larger 
grains and outward transport of small grains lead to the fact that 
the dust-to-gas ratio is reduced by 0.5-1.5 orders of magnitude 
compared to a stationary model. This effect is also discussed in 
Section 3.4. 



3.3. Fiducial model with fragmentation 

The situation changes significantly, if we take grain fragmen- 
tation into account. As discussed in Section 2.6.2, we consider 
two different kinds of outcomes for grain-grain collisions with 
relative velocities larger than the fragmentation velocity uf. cra- 
tering (if the masses differ by less than one order of magnitude) 
and complete fragmentation (otherwise). 

For particles larger than about St ss 10~ 3 , relative velocities 
are dominated by turbulent motion (and to a lesser extend by ver- 
tical settling). Since the relative velocities increase with Stokes 
number (and therefore with grain size), we can calculate the ap- 
proximate position of the fragmentation barrier by equating the 
assumed fragmentation velocity Uf with the approximate relative 



10 T. Bimstiel et al.: Gas- and dust evolution in protoplanetary disks 



1 2 1 2 

10 10 10 10 10 10 




o 1 2 1 2 

10 10 10 10 10 10 

r [AU] r [ATI] 



Fig. 5. Snapshots of the vertically integrated dust density distributions (defined in Eq. 37) of a steady state disk as in B08a (left 
column) and of an evolving disk (fiducial model, right column). No coagulation is calculated within the evaporation radius (denoted 
by the dash-dotted line), fragmentation is not taken into account in both simulations. The solid line shows the particle size corre- 
sponding to a Stokes number of unity. Since «st=i K 2 g (see Eq. 14) this curve in fact has the same shape as E g (r), so it reflects, as 
a "bonus", what the gas disk looks like. The radius dividing the evolving disk into accreting and expanding regions is marked by 
the dotted line and the arrows. Particles which are located below the dashed line have a positive flux in the radial direction due to 
coupling to the expanding gas disk and turbulent mixing (particles within the closed contour in the upper right plot have an inward 
pointing flux). 
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Fig. 6. Comparison of the radial dependence of the dust-to-gas ratio for the stationary simulations (thick lines) and the evolving disk 
simulations (thin lines). The four panels compare the results after 1 Myr of evolution with/without fragmentation (left/right column) 
and with/without effects of non-axisymmetric instabilities (top/bottom row). It can be seen that the dust-to-gas ratio of the evolving 
disk simulations is almost everywhere lower than the one of the stationary simulations. See Section 3.4 for more details. 



velocities of the particles, 



Particles larger than this size are subject to high-velocity colli- 
sions which will erode or completely fragment those particles. 
This is only a rough estimate as the relative velocities also de- 
pend on the size of the smaller particle and radial drift also influ- 
ence the grain distribution. However Eq. 52 reproduces well the 
upper end of the size distribution in most of our simulations and 
therefore helps understanding the influence of various parame- 
ters on the outcome of these simulations. 

The evolution of the grain size distribution including frag- 
mentation is depicted in Figure 7. The initial condition is quickly 
forgotten since particles grow on very short timescales to sizes at 
which they start to fragment. The resulting fragments contribute 
again to the growth process until a semi-equilibrium of growth 
and fragmentation is reached. 

It can be seen that particles stay much smaller than in 
the model without fragmentation. This means that they are 
less affected by radial drift on the one hand and better trans- 
ported along with the expanding gas disk on the other hand. 
Consequently, considerable amounts of dust can reach radii of 
the order of 100 AU, as seen in Figure 8. 



The approximate maximum Stokes number, defined in 
Eq. 52, is inversely proportional to the temperature (since rel- 
ative velocities are proportional to c s ), which means that in re- 
gions with lower temperature, particles can reach larger Stokes 
numbers. By equating drift and drag velocities of the particles 
(cf. Eq. 19), it can be shown that the radial velocities of particles 
with Stokes numbers larger than about a/2, are being dominated 
by radial drift. 

Due to the high temperatures below ~5 AU (caused by vis- 
cous heating), St max stays below this value which prevents any 
significant radial drift within this radius, particles inside 5 AU 
are therefore only transported along with the accreting gas. 
Particles at larger radii and lower temperatures can drift (al- 
though only slightly), which means that there is a continuous 
transport of dust from the outer regions into the inner regions. 
Once these particles arrive in the hot region, they get "stuck" 
because their Stokes number drops below a/2. The gas within 
about 5 AU is therefore enriched in dust, as seen in Figure 8. 
The dust-to-gas ratio at 1 AU after 1 Myr is increased by 25% 
over the value of in-falling matter, which is taken to be 0.01. 

Figure 8 also shows a relatively sharp decrease in the dust to 
gas ratio at a few hundred AU. At these radii, the gas densities 
become so small that even the smallest dust particles decouple 
from the gas and start to drift inwards. 
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The thick line in Figure 8 shows as comparison the dust- 
to-gas ratio of the stationary disk model (cf. left column of 
Figure 5) after 1 Myr, which starts with a radially constant initial 
dust-to-gas ratio of 0.01. 

Figure 9 shows the semi-equilibrium grain surface density 
distribution at 1, 10 and 100 AU in the fiducial model with frag- 
mentation at 1 Myr. The exact shape of these distributions de- 
pends very much on the prescription of fragmentation and cra- 
tering. In general the overall shape of these semi-equilibrium 
distributions is always the same: a power-law or nearly constant 
distribution (in <x) for small grains and a peak at some grain size 
flmax, beyond which the distribution drops dramatically. The peak 
near the upper end of the distribution is caused by cratering. This 
can be understood by looking at the collision velocities: the rel- 
ative velocity of two particles increases with the grain size but 
it is lower for equal-sized collisions than for collisions with par- 
ticles of very different sizes (see Figure 2). The largest parti- 
cles in the distribution have relative velocities with similar sized 
particles which lie just below the fragmentation velocity (other- 
wise they would fragment). This means that the relative veloci- 
ties with much smaller particles (which are too small to fragment 
the bigger particles but can still damage them via cratering) are 
above this critical velocity. This inhibits the further growth of 
the big particles beyond a max , causing a "traffic jam" close to 
the fragmentation barrier. The peak in the distribution represents 
this traffic jam. 
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Fig. 8. Evolution of the radial dependence of the dust-to-gas ra- 
tio in the fiducial model including fragmentation with the times 
corresponding to the snapshots shown in Figure 7. The initial 
dust-to-gas ratio is taken to be 0.01. The thick dashed curve rep- 
resents the result at 1 Myr of the static disk model for compari- 
son. 



Fig. 7. Evolution of the dust density distribution of the fidu- 
cial model as in Figure 5 but with fragmentation included. The 
dashed contour line (in the lower two panels) around the upper 
end of the size distribution and around small particles at > 60 AU 
marks those parts of the distribution which have a positive (out- 
ward pointing) fluxes. 



3.4. Influences of the infall model 

In the fiducial model without fragmentation, continuous resup- 
ply of material by infall is the cause why the disk has much 
more small grains than compared to the stationary disk model 
(cf. Figure 5), which relatively quickly consumes all available 
micrometer sized dust. The effect has already been found in 
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Fig. 9. Vertically integrated (cf. Eq. 37) grain surface density 
distributions as function of grain radius at a distance of 1 AU 
(solid), 10 AU (dashed) and 100 AU (dot-dashed) from the 
star. These curves represent slices through the bottom panel of 
Figure 7. 



Dominik & Dullemond (2008): if all grains start to grow at the 
same time, then the bulk of the mass grows in a relatively thin 
peak to larger sizes (see Figure 6 in B08a). However if the bulk 
of the mass already resides in particles of larger size, then addi- 
tional supply of small grains by infall is not swept up effectively 
because of the following: firstly, the number density of large par- 
ticles is small (they may be dominating the mass, but not neces- 
sarily the number density distribution) and secondly, they only 
reside in a thin mid-plane layer while the scale height of small 
particles equals the gas scale height. 

We studied how much the disk evolution depends on changes 
in the infall model. 

For a given cloud mass, the so called centrifugal radius r centr , 
which was defined in Eq. 28, depends on the temperature and 
the angular velocity of the cloud. Both can be varied result- 
ing in a large range of possible centrifugal radii reaching from 
a few to several hundred AU. Since the centrifugal radius is 
the relevant parameter, we varied only the rotation rate of the 
cloud core. We performed simulations with three different rota- 
tion rates which correspond to centrifugal radii of about 8 (fidu- 
cial model), 33 AU, and 100 AU. For each centrifugal radius, we 
performed two simulations: one which includes effects of gravi- 
tational instabilities (GI) - i.e. increased a during infall and ac- 
cording to Eq. 6 - and one which neglects them. 

However for a centrifugal radius of 100 AU, too much mat- 
ter is loaded onto the cold outer parts of the disk and conse- 
quently, the disk would fragment through gravitational instabil- 
ity. We cannot treat this in our simulations, hence, for the case 
of 100 AU, we show only results which neglect all GI effects. 

The resulting dust-to-gas ratios are being shown in Figure 6. 

Two general aspects change in the case of higher rotation 
rates: firstly, more of the initial cloud mass has to be accreted 
onto the star by going through the outer parts of the disk. 
Consequently, the disk is more extended and more massive than 
compared to the case of low rotation rate. 

Secondly, as aforementioned the high surface densities in the 
colder regions at larger radii cause the disk to become less grav- 
itationally stable. 



If grain fragmentation is not taken into account in the sim- 
ulations, both effects cause more dust mass to be transported to 
larger radii. Growth and drift timescales are increasing with ra- 
dius and the dust disk with centrifugal radius of 33 AU (100 AU) 
can stay 5 (35) times more massive than in the low angular mo- 
mentum case after 1 Myr if GI effects are neglected. 

If GI effects are included, matter is even more effectively 
transported outward, the dust-to-disk mass ratio for 8 and 33 AU 
is increased from 5 to 8. 

However if fragmentation is included, it does not matter so 
significantly, where the dust mass is deposited onto the disk 
since grains stay so small during the build-up phase of the disk 
(due to grain fragmentation by turbulent velocities) that they are 
well coupled to the gas. Outwards of ~ 10 AU (without GI ef- 
fects) or of a few hundred AU (if GI effects are included), the 
gas densities become so small that even the smallest grains start 
do decouple from the gas. They are therefore not as effectively 
transported outwards. In these regions, the amount of dust de- 
pends on the final centrifugal radius while at smaller radii, tur- 
bulent mixing quickly evens out all differences in the dust-to-gas 
ratio (see left column of Figure 6). 

It can be seen, that in all simulations, the dust-to-gas ratio is 
lower than in the stationary disk model. The trend in the upper 
right panel in Figure 6 suggests that for a centrifugal radius of 
100 AU and the enhanced radial transport by GI effects might 
have a higher dust-to-gas ratio than the stationary disk model. 
However in this case, the disk would become extremely unstable 
and would therefore fragment. 

The reason for this is the following: to be able to compare 
both simulations, the total mass of the disk-star system is the 
same as in the stationary disk models. How the total mass is 
distributed onto disk and star depends on the prescription of 
infall. If a centrifugal radius of 100 AU is used, the disk be- 
comes so massive that it significantly exceeds the stability crite- 
rion AW/M* < 0.1. 

3.5. The radial drift barrier revisited 

According to the current understanding of planet formation, sev- 
eral mechanisms seem to prevent the formation of large bodies 
via coagulation quite rigorously. The most famous ones - radial 
drift and fragmentation - have already been introduced above. 
Radial drift has first been discussed by Weidenschilling (1977), 
while the importance of the fragmentation barrier (which may 
prevent grain growth at even smaller sizes) was discussed in 
B08a. In the following, we want to test some ideas about how to 
weaken or to overcome these barriers apart from those already 
studied in B08a. 

B08a has quantified the radial drift barrier by equating the 
timescales of growth and radial drift which leads to the condition 



GO 



1 

r' 



(53) 



where eo is the dust-to-gas ratio and and r g are the drift and 
growth timescales respectively. The parameter y describes how 
much more efficient growth around St=l must be, so that the 
particles are not removed by radial drift. To overcome the drift 
barrier, obviously either particle growth must be accelerated, or 
the drift efficiency has to be decreased. B08a have numerically 
measured the parameter y to be around 12. In other words, this 
means that the growth timescales have to be decreased (e.g. by 
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Fig. 10. Evolution of the dust surface density distribution of the 
fiducial model at 200000 years for different drift efficiencies E^, 
without fragmentation. The solid line denotes the grain size ast=i 
of particles with a Stokes number of unity. Gas outside of the 
radius denoted by the dotted line as well as particles below the 
dashed line have positive radial velocities. See section 3.2 for 
more discussion. 



an increased dust-to-gas ratio) until the condition in Eq. 53 is 
fulfilled. 

However, there are other ways of breaking through the drift 
barrier. Firstly, the drift timescale for St=l particles also depends 
on the temperature (via the pressure gradient). A simple approx- 



imation from Eq. 53 with a 0.5 M star and a dust-to-gas ratio of 
0.01 gives 

r<103K(^)"\ (54) 

which means that particles should be able to break through the 
drift barrier at 1 AU if the temperature is below ~100 K (or 
200 K for a solar mass star). Dullemond et al. (2002) have con- 
structed vertical structure models of passively irradiated circum- 
stellar disks using full frequency- and angle-dependent radiative 
transfer. They show that the mid-plane temperature of such a T 
Tauri like system at 1 AU can be as low as 60 K. Reducing the 
temperature by some factor reduces the drift time scale by the a 
factor of similar size which we will call the radial drift efficiency 
Ei (cf. Eq. 20). 

Zonal flows as presented in Johansen et al. (2006) could be 
an alternative way of decreasing the efficiency of radial drift av- 
eraged over typical time scales of particle growth. Johansen et al. 
(2006) found a reduction of the radial drift velocity of up to 40% 
for meter-sized particles. 

Meridional flows (e.g., Urpin 1984; Kley & Lin 1992) might 
also seem interesting in this context, however they do not di- 
rectly influence the radial drift efficiency but rather reverse the 
gas-drag effect. This might be important for small particles 
(which, however are not strongly settling to the mid-plane) but 
for St=l particles, a would have to be extremely high to have 
significant influence: even a = 0. 1 would result in a reduction of 
the particles radial velocity by approximately only a few percent. 

Another possibility to weaken the drift barrier is changing its 
radial dependence. The reason for this is the following: particle 
radial drift is only a barrier if it prevents particles to cross the 
size flst=i- Since particles grow while drifting, the particle size 
corresponding to St= 1 needs to increase as well, to be a barrier. 
Otherwise drifting particles would grow (at least partly) through 
the barrier while they are drifting. If as t =i is decreasing in the 
direction toward the star, then a particle that drifts inwards would 
have an increasing Stokes number even if the particle does not 
grow at all. 

In the Epstein regime, the size corresponding to St= 1 is pro- 
portional to the gas surface density 



ast=i = 



22g 



(55) 



meaning that a relatively flat profile of surface density (or even a 
profile with positive slope) is needed to allow particles to grow 
through the barrier. However, our simulations of the viscous gas 
disk evolution does not yield surface density profiles with posi- 
tive slopes outside the dust evaporation radius. 

To quantify the arguments above, we have performed sim- 
ulations with varying drift efficiency to test how much the 
radial drift has to be reduced to allow break through. We have 
additionally included the first Stokes drag regime to see how the 
radial drift of particles is influenced by it. 

Figure 10 shows the grain surface density distribution after 
200000 years of evolution for three different drag efficiencies. 
The most obvious changes can be seen in the region where the 
flst=i line (solid line) is relatively flat: the grain distribution is 
shifted towards larger Stokes numbers. As explained above, par- 
ticles can grow while drifting, which can be seen in the case 
of E,\ = 0.5. The Stokes number and size of the largest parti- 
cles is significantly increasing towards smaller radii. However 
the radial drift efficiency has to be reduced by 80% to produce 
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Fig. 12. Dust surface density distributions (top row) and the according solid-to-gas ratio (bottom row) for the case of radius- 
dependent fragmentation velocity after 2xl0 5 years of evolution. In the upper row, the vertical dashed line denotes the position 
of the snow line at the mid-plane, the solid line corresponds to as t =i and the dot-dashed line shows the approximate location of 
the fragmentation barrier according to Eq. 52. The snow line on the right plot lies further outside since viscous heating is stronger 
if a is larger. In the bottom row, the solid line denotes the vertically integrated dust-to-gas ratio while the dashed line denotes the 
dust-to-gas ratio at the disk mid-plane. The icy dust grains outside the snow line are assumed to fragment at a critical velocity of 
10 ms"' while particles inside the snow line fragment at 1 m s . The plots on the left and right hand side differ in the amount of 
turbulence in the disk (a = 10~ 3 and 10~ 2 , respectively). 



particles which are large enough to escape the drift regime and 
are therefore not lost to the star. 

The situation changes, if the Stokes drag is taken into ac- 
count: if gas surface densities are high enough for the dust par- 
ticles to get into a different drag regime, a change in the radial 
dependency of ast=i can be achieved. The Epstein drag regime 
for particles with Stokes number of unity is only valid if 



£ g < 108- 



\200K7 VAU/ 



to*) 


4 











(56) 



otherwise, drag forces have to be calculated according to the 
Stokes drag law since the Knudsen number becomes smaller 
than 4/9 (see Weidenschilling 1977). The Stokes number is then 
given by 



St = 



V2^ 



Ps cr H , a* 



h: 



H nip 



(57) 



with <th, being the cross section of molecular hydrogen. 
Interestingly, ast=i is independent of Z g and proportional to the 
square root of the pressure scale height which decreases towards 
smaller radii. This means that - as long as the surface density 
is high enough - it does not depend on the radial profile of the 
surface density. In this regime the radial drift itself could move 
particles over the drift barrier since drifting inwards increases 
the Stokes number of a particle without increasing its size. 

Results of simulations which include the Stokes drag are 
shown in Figure 11. In this case, particles can already break 
through the drift barrier if E4 < 0.75. This value and the posi- 
tion of the breakthrough depends on where the drag law changes 
from Epstein to Stokes regime and therefore on the disk surface 
density. As noted above, larger surface densities generally shift 
the position of regime change towards larger radii making it eas- 
ier for particles to break through the drift barrier. 
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Fig. 11. Dust grain surface density distribution as in Figure 10 
at 200 000 years but including the Stokes drag regime. The drift 
efficiency is set to = 0.75 and fragmentation is not taken 
into account. It can be seen that «st=i (solid line) increases with 
radius until about 1 AU, which facilitates the break through the 
drift barrier. 

It should be noted that the physical way to avoid the Stokes 
drag regime is to decrease the surface densities, however we 
chose the same initial conditions for both cases - with and with- 
out Stokes drag - and just neglected the Stokes drag in the latter 
computations to be able to compare the efficiency factors inde- 
pendent of other parameters such as disk mass or temperature. 

3.6. The fragmentation barrier revisited 

In the previous section, we have shown that several mechanisms 
could allow particles to break through the radial drift barrier, 
however fragmentation puts even stronger constraints on the for- 
mation of planetesimals. 

As shown by Ormel & Cuzzi (2007), the largest relative ve- 
locities are of the order of 

Aw max - V2ac s . (58) 

If particles should be able to break through the fragmentation 
barrier, then they need to survive these large relative velocities, 
meaning that AM max has to be smaller than the fragmentation ve- 
locity of the particles, or 

— > V2^. (59) 

c s 

This condition is hard to fulfill with reasonable fragmentation 
velocities, unless a is very small. E.g., for a = 10~ 5 and a tem- 
perature of 200 to 250 K, the fragmentation velocity needs to 
be higher than 4 m s~', which could already seen in the sim- 
ulations by Brauer et al. (2008b), who have simulated particle 
growth near the snow-line. 

Even in the case of very low turbulence, relative azimuthal 
velocities of large (St > 1) and small grains (St «: 1) are of the 
order of 30 m s -1 , which means that large particles are constantly 
being 'sand-blasted' by small grains. The only way of reducing 
these velocities significantly is decreasing the pressure gradient 
(see Equations 20 and 46). 



Another possibility to overcome this problem would be 
if high-velocity impacts of smaller particles would cause net 
growth, as has been found experimentally by Wurm et al. (2005) 
and Teiser & Wurm (2009). 

Taken together, these facts make environments as the in- 
ner edge of dead zones ideal places for grain growth (see 
Brauer et al. 2008b; Kretke & Lin 2007): shutting of MRI leads 
to low values of a, which are needed to reduce turbulent relative 
velocities and the low pressure gradients prevent radial drift and 
high azimuthal relative velocities. 

3. 7. Dust enhancement inside the snow line 

As already noted in a previous paper (Birnstiel et al. 2009), sig- 
nificant loss of dust by radial drift can be prevented by assuring 
that particles stay small enough and are therefore not influenced 
by radial drift. For typical values of a (10~ 3 - 10~ 2 ), this means 
that the fragmentation velocity must be smaller than about 0.5- 
5 ms 4 . If particles have higher tensile strength, they can grow 
to larger sizes which are again affected by radial drift. 

Typical fragmentation velocities for silicate grains deter- 
mined both theoretically and experimentally are of the order of a 
few m s~' (for a review, see Blum & Wurm 2008). The composi- 
tion of particles outside the snow-line is expected to change due 
to the presence of ices. This can influence material properties 
and increase the fragmentation velocity (see Schafer et al. 2007; 
Wada et al. 2009). 

We have performed simulations with a radially varying frag- 
mentation speed. We assume the fragmentation speed inside the 
snow-line to be 1 m s , outside the snow-line tobe 10 ms -1 . 
It should be noted that we do not follow the abundance of water 
in the disk or the composition of grains, we only assume parti- 
cles outside the snow line to have larger tensile strength due to 
the presence of ice. To be able to compare both simulations, we 
used the same 1 Myr old 0.09 M Q gas disk around a solar mass 
star as initial condition. The gas surface density profile of this 
disk was derived by a separate run of the disk evolution code. 
We used this gas surface density profile and a radially constant 
dust-to-gas ratio as initial condition for the simulations presented 
in this subsection. Apart from the level of turbulence, the input 
for both simulations is identical, the results are therefore com- 
pletely independent of uncertainties caused by the choice of the 
infall model. 

Results of the simulations are shown in Figure 12. A one or- 
der of magnitude higher fragmentation velocity causes the max- 
imum grain size to be about two orders of magnitude larger, 
which follows from Eq. 52 since St max oc a max in the Epstein 
regime (all particles in these simulations are small enough to 
be in the Epstein regime). This effect can be seen in Figure 12. 
Particles outside the snow-line are therefore more strongly drift- 
ing inwards (because they reach larger Stokes numbers) where 
they are being pulverized as soon as they enter the region within 
the snow-line. 

Strong drift outside the snow line and weaker radial drift in- 
side the snow line cause the dust-to-gas ratio within the snow 
line to increase significantly (see bottom row of Figure 12): in 
the case of a — 10 -3 , the dust-to-gas ratio reaches values be- 
tween 0.39 and 0. 10 in the region from 0.2 to 1 .9 AU. 

Simulations for a less massive star (0.5 M ) show the same 
behavior, however the increase in dust-to-gas ratio is not as high 
as for a solar mass star (dust-to-gas ratio of 0.27-0.20 from 0.2- 
4 AU). 

This effect is not as significant in the case of stronger tur- 
bulence, where the maximum dust-to-gas ratio is around 0.027. 
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The evolution of the dust-to-gas ratio at a distance of 1 AU from 
the star is plotted for both cases in Figure 13 (the minor bump is 
an artifact of the initial condition). 

The reason for this difference lies in the locations of the 
drift and fragmentation barriers. The approximate position of 
the fragmentation barrier is represented by the dot-dashed line in 
Figure 12. The radial drift barrier cannot be defined as sharply, 
however radial drift is strongest at a Stokes number of unity, 
which corresponds to the solid line in Figure 12. An increase of 
a by one order of magnitude lowers the fragmentation barrier by 
about one order of magnitude in grain size. 

In the lower turbulence case, the fragmentation barrier lies 
close to flst=i- Most particles are therefore drifting inwards be- 
fore they are large enough to experience the fragmentation bar- 
rier. Hence, the outer parts of the disk are significantly depleted 
in small grains. 

In contrast to this case, fragmentation is the stronger barrier 
for grain growth throughout the disk in the high turbulence simu- 
lation (right column in Figure 12). It can be seen that particles of 
smaller sizes are constantly being replenished by fragmentation. 

With these results in mind, the evolution of the disk mass 
(bottom panel of Figure 14) seems counter-intuitive: the mass 
of the high turbulence dust disk is decaying faster than in the 
low turbulence case. This can be understood by looking at the 
global dust-to-gas ratio of the disks (top panel of Figure 14) 
which does not differ much in both cases. This means that the 
increased dust mass loss in the high turbulence disk is due to the 
underlying evolution of the gas disk. Particles in the high turbu- 
lence disk may have smaller Stokes numbers (causing drift to be 
less efficient), however the inward dragging by the accreting gas 
is stronger in this case. 

To show how much the dust evolution depends on the frag- 
mentation velocity, we included the case of a lower fragmenta- 
tion velocity throughout the disk in Figure 14. It can be seen 
that the dust mass is retained at its initial value for much longer 
timescales. 




time [yr] 



Fig. 13. Dust-to-gas ratio at a distance of 1 AU from the central 
star as a function of time for the case of low (a = 10~ 3 , solid 
line) and high (a = 1CT 2 , dashed line) turbulence. 
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Fig. 14. Evolution of the global dust-to-gas ratio (top panel) and 
the dust disk mass (bottom panel) for the simulations shown in 
Figure 12. The solid and dashed lines correspond to the low and 
high turbulence case, respectively. The dash-dotted line shows 
the evolution of the disk if a low fragmentation velocity is as- 
sumed throughout the disk. 

4. Discussion and conclusions 

We constructed a new model for growth and fragmentation of 
dust in circumstellar disks. We combined the (size and radial) 
evolution of dust of B08a with a viscous gas evolution code 
which takes into account the spreading and accretion, irradiation 
and viscous heating of the gas disk. The dust model includes the 
growth/fragmentation, radial drift/drag and radial mixing of the 
dust. We re-implemented and substantially improved the numer- 
ical treatment of the Smoluchowski equation of B08a to solve for 
the combined size and radial evolution of dust in a fully implicit, 
un-split scheme (see Appendix A). In addition to that, we also in- 
cluded more physics such as relative azimuthal velocities, radial 
dependence of fragmentation critical velocities and the Stokes 
drag regime. The code has been tested extensively and was found 
to be very accurate and mass-conserving (see Appendix B). 

We compared our results of grain growth in evolving proto- 
planetary disks to those of steady state disk simulations, similar 
to B08a. In spite of many differences in details, we confirm the 
most general result of B08a: radial drift and particle fragmenta- 
tion set strong barriers to particle growth. If fragmentation is in- 
cluded in the calculations, then it poses the strongest obstacle for 
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the formation of planetesimals. Very low turbulence (a < 1CT 5 ) 
and fragmentation velocities of more than a few m s are needed 
to be able to overcome the fragmentation barrier in the case of 
turbulent relative velocities. 

This model includes also the initial build-up phase of the 
disk, which is still a very poorly understood phase of disk evo- 
lution. We use the Shu-Ulrich infall model which represents a 
strong simplification. However, the following novel findings of 
this work do not or only weakly depend on the build-up phase of 
the disk: 

- Apart from an increased dust-to-gas ratio (as in B08a), other 
mechanisms such as streaming instabilities or a decreased 
temperature may be able to weaken this barrier by decreas- 
ing the efficiency of radial drift. We found that in simula- 
tions without fragmentation the radial drift efficiency needs 
to be reduced by 80% to produce particles which crossed the 
meter-size barrier and are large enough to resist radial drift. 

- If the gas surface density is above a certain threshold (in our 
simulations about 140 g ctrT 2 at 1 AU or 330 g crrT 2 at 5 AU, 
see Eq. 56) then the drag force which acts on the dust par- 
ticles has to be calculated according to the Stokes drag law, 
instead of the Epstein drag law. The drift barrier in this drag 
regime is shifting to smaller sizes for smaller radii (inde- 
pendent on the radial profile of the surface density) which 
means that pure radial drift can already transport dust grains 
over the drift barrier or at least to larger Stokes numbers even 
without simultaneous grain growth. In this case, the drift ef- 
ficiency needs to be reduced only by about 25% to produce 
large bodies. 

- If relative azimuthal velocities are included, then grains with 
St > 1 are constantly 'sand-blasted' by small grains (if 
they are present) which (in our prescription of fragmen- 
tation) causes erosion and stops grain-growth even in the 
case of low turbulence. Only decreasing the radial pres- 
sure gradient significantly weakens both relative azimuthal 
and radial velocities. Low turbulence and a small radial 
pressure gradient together are needed to allow larger bod- 
ies to form. These conditions may be fulfilled at the inner 
edge of dead zones (Brauer et al. 2008b; Kretke & Lin 2007; 
see also Dzyurkevich et al., in press). Future work needs 
to investigate the disk evolution and grain growth of disks 
with dead zones. Our prescription of fragmentation and ero- 
sion may also need rethinking since Wurm et al. (2005) and 
Teiser & Wurm (2009) find net-growth by high velocity im- 
pacts of small particles onto larger bodies. 

- Higher tensile strengths of particles outside the snow-line 
allows particles to grow to larger sizes, which are more 
strongly affected by radial drift. Particles therefore drift from 
outside the snow-line to smaller radii where they fragment 
and almost stop drifting (since their radial velocity is de- 
creased by almost two orders of magnitude). This can cause 
an increase the dust-to-gas ratio inside the snow-line by more 
than 1.5 orders of magnitude. 

- The critical fragmentation velocity and its radial dependence 
(and to a lesser extent the level of turbulence) is a very im- 
portant parameter determining if the dust disk is drift or frag- 
mentation dominated. A drift dominated disk is significantly 
depleted in small grains and only a fragmentation dominated 
disk can retain a significant amount of dust for millions of 
years as is observed in T Tauri disks. 

The following results depend on the build-up phase of the 
disk. However unless the collapse of the parent cloud is not 



inside-out or so fast to cause disk fragmentation, we expect only 
slight alteration of the results: 

- Disk spreading causes small particles (< 10 fim) to be trans- 
ported outward at radii of ~60-190 AU even in 1 Myr old 
disks. 

- Small particles provided by infalling material are not effec- 
tively swept up by large grains if the bulk of the dust mass 
has already grown to larger sizes. 

- In an inside-out build-up of circumstellar disks, grains 
growth is very fast (timescales of some 100 years) because 
densities are high and orbital timescales are small. Large 
grains are quickly lost due to drift towards the star if frag- 
mentation is neglected. Fragmentation is firstly needed to 
keep grains small enough to be able to transport a signif- 
icant amount of dust to large radii by disk spreading and 
secondly to retain it in the disk by preventing strong radial 
inward drift. 
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Working out these equations and separating the values of N 
leads to 

K = a„ ■ a£1 + b„ ■ n;; 1 + c„ ■ + d„ (a.8) 



with the coefficients 



max(0, u n _i) ■ S n _i + 



B„ = i-AtL n +mmax(0,u n+k )-S, 



-min(0, m„_i) ■ S n _i. 

K±i s» D ,\„ i ,/! „ i ' s „ >*> 



(A„ +1 -*„)/!„ 



(X„-A„_,)/l„ 



C„ = ^|min(0, M „ + .)-5„ 



(i„ + i-Jc„)A„ +1 



D n = -At-K„. 



(A.9) 



Appendix A: Algorithm 

In the next two sections, we will first discuss how the equations 
of radial evolution of gas (Eq. 1) and dust (Eq. 21) as well as the 
coagulation/fragmentation of dust (Eq. 35) are solved separately. 
In Section A. 3 we will then describe how this model treats the 
radial and the size evolution of dust in an un-split, fully implicit 
way. 



A. 1 . Advection-diffusion Algorithm 

To be able to also model both, the evolution of dust and gas 
implicitly, we constructed a scheme which solves a general form 
of an advection-diffusion equation, 



Eq. A.8 can now be solved by any matrix-solver, but since it is 
a tri-diagonal matrix, the fastest analytical way to solve it is by 
forward elimination/backward substitution. 

It should be noted that Eq. A. 1 is implicit only in TV which 
means that the equations we solve are only implicit in the sur- 
face density. In the case of the viscous accretion disk, described 
by Eq. 1, we face the problem that also the turbulent gas viscos- 
ity v g depends on the temperature which (in the case of viscous 
heating) depends on the surface density. This can cause numeri- 
cal instabilities to develop. 

To stabilize the code, we use a scheme which estimates the 
temperature in several predictor steps. The actual time step is 
then done with the predicted temperature. 



dN 



d_ 
dx 



(N-u)- 



d 








dx 






= K + L-N (A.l) 



which can be adapted to both, Eq. 1 and Eq. 21 by proper choice 
of parameters u, D^, g, h, K and L. 

We use a flux-conserving donor-cell scheme which is im- 
plicit in N. The time derivative in Eq. A.l, written in a dis- 
cretized way becomes 



dt t i+ i - U 



(A.2) 



where ; denotes time-dimension and n denotes space-dimension. 
The advective part is discretized as 



dx 

where F' + \ and S 



(N-u) 



pi+l 

»+4 



(A.3) 



V„ V n 
are the future flux and the surface between 



cell n and cell n + 1 and V„ is the volume of cell «. The advective 
interface fluxes can then be written as 



N„ ■ max(0, u n+ \_ ) + N„+i ■ min(0, u n+ \_ ) 



F' + \ = N„-i ■ max(0, u n _\_ ) + 7V„ ■ min(0, u n _\_ ) 



The diffusive interface flux becomes 



7 i+ 1 
d,n- 



= D A 



(A.4) 
(A.5) 

(A.6) 
(A.7) 



A.2. Coagulation Algorithm 

Discretizing Eq. 35 on a mass grid m,- gives 

a 



jn k (r, z) = 2 M ijk rii(r, z) n/r, z), (A. 10) 



dt 

•j 

where the dust particle number density is 



«;(>, z) 



n 



(m, r, z) dm 



(A. 11) 



If we assume that the coagulation and fragmentation kernels 
are constant in z and that the vertical distribution of grains is 
a Gaussian with a scale height according to Eq. 51, 



nk(r,z) = 



N k (r) 



exp - 



2h k (ry 



(A. 12) 



we can vertically integrate the coagulation/fragmentation equa- 
tion. 



at Nk{r) = L hkdz 
u 

( h}+h 2 : 



hi hj 



exp 



-^^■z 2 \dz 

2h]h) 



(A. 13) 



YjMijkNMNjir), 
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where 



Mi 



ijk 



M ijk . 



(A. 14) 



Discretizing Eq. A. 13 also in radial direction and rewriting it in 
terms of the quantity 



yields 



M/ = N k (r,) ■ r, 



d v~i M ijk 

~Q t N u = 2~^r ' Nil ' Njl := Su ' 

ij 



(A. 15) 



(A. 16) 



where the vector S = {Ski}k=\,m is the source function for each of 
the m mass bins. 

The numerical change of Nu within a time step Af = f, - f,_i 
is AW/;./ = A/It 1 - N' k[ . The time-discretized version of Eq. A. 16 
then becomes (omitting second order terms) 

u 



=Tj~t ( NaNji + ANaNji + NuANji + ^ NhAN ^) ■ 

(A. 17) 



Since the first term on the right hand side is the explicit source 
function, we can write 



AM/ 
Af 



M ijk + M jik 

n 



Nji ■ AM/- (A. 18) 



Using the vectors AN = {AJV}/t=i,n m and N = {N]k=\,n m , this can 
be rewritten in matrix notation, 



ih-') 



AN = S. 



Were 



Mjjk + Mj^ 

n 



N 



(A.19) 



(A.20) 



denotes the Jacobian of the source function and 1 the unity ma- 
trix. The solution for the future values can now be derived by 
inverting the matrix in Eq. A. 19, 



N'' +1 = N' + AN 



MlH" 



(A.21) 



A3. Fully implicit scheme for radial motion and coagulation 

To be able to solve the radial motion and the Smoluchowski 
equation fully implicitly, we rewrite Equation A. 8 as 



M • N i+1 = N - D, 



(A.22) 



where M is the tri-diagonal matrix with entries A, B and C and 
source term D which are given by Eq. A. 9. M is now rewritten 
by separating off the diagonal terms and by dividing by Af, 



Fig.A.l. Pictographic representation of the matrix on the left 
hand side of Eq. A. 26 with six radial and five mass grid points. 
White elements are zero, dark grey elements contain contribu- 
tions from radial transport of dust (J), light gray elements con- 
tain contributions from the coagulation/fragmentation (M), and 
black elements contain both contributions. The upper left and 
lower 5 rows represent the boundary conditions where coagula- 
tion/fragmentation is not taken into account. The matrix in the 
simulations would typically have a size of 15000 2 . 



which brings Eq. A.22 in a form similar to Eq. A.19, 



(£♦*)■ 



AN 



-M ■ N - D/Af 



(A.24) 



The coagulation/fragmentation is now determined by the matrix 
J and the source vector S and similarly, the radial motion is de- 
termined by matrix M and source vector 



R = -M ■ N - D/Af. 



(A.25) 



This allows us to combine both operators into one matrix equa- 
tion, 



+ M - jj ■ AN = R + S. 



(A.26) 



In principle, to solve for the vector AN, the matrix on the left 
hand side of Eq. A.26 has to be inverted. This is a numeri- 
cally challenging task since the inverse matrix in our simulations 
would have about 150-500 million elements. The matrix on the 
left hand side of Eq. A.26, however is a sparse matrix (schemat- 
ically depicted in Figure A.l). 

We can therefore solve Eq. A.26 by an iterative incomplete 
LU decomposition solver for sparse matrices provided by the 
Sixpack library 2 of S. E. Norris. 



Appendix B: Test cases 

To check if the numerical implementation presented above ac- 
curately solves Eq. A.26, we compare results of the simulation 
to some analytical solutions: The growth rates of particles can 
be approximated if we assume the grain size distribution to be 
a delta function and the sticking probability to be unity. The in- 
crease of mass per collision is then given by the mass of the 



M = Af ■ M + 1 



(A.23) 



available from www . engineers . auckland . ac . nz/~ snor007 
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Fig. B.l. Test case for only radial drift, without coagulation. The 
solid line denotes the mass-averaged position of the radial distri- 
bution of grains for each grain size after 10 3 years. The dashed 
line is the expected solution for a single particle, the dotted line 
denotes the initial radius of the particles. 
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Fig.A.2. Comparison between simulation result and analytical 
solution for growth of equal sized particles. The solid lines de- 
note the position of the peak of the grain size distribution. In 
the top panel, only Brownian motion is considered as source of 
relative particle velocities, in the bottom panel, turbulent rela- 
tive velocities are considered as well. The parameters of these 
simulation are T = 196 K, p s =1.6 g/cnr, S a = 18 g/cm 2 , 
S d = 0.18 g/cm 2 and a = 1(T 3 . 



particle, m divided by the time between two collisions, r, which 
can be written as 



dm m 
dt t ' 

Using dm = Anp^da and r = m/(pd cr Am), we derive 

? = ^ Ah. 
dt p s 



(B.l) 



(B.2) 



where Am is the relative velocity, pa the dust density, p s the solid 
density of the dust particles and cr = Ana 2 the cross section of 
the collision. 

With this formula, we can estimate the growth time scales if 
the relative velocities of the dust grains is given. For equal sized 
particles and Brownian motion, we get 



Am 

for turbulent motion, 
Ormel & Cuzzi 2007) 



BM 



16 k h T 



(B.3) 



the relative velocities are (see 



Am 



Integrating Eq. B.2 give 
a(t) 



( c s V2 a St for St «: 1 
H forSt»l 




A (f-f ) + fl n 2 



for Brownian motion and 



a(t) 



a ■ exp ^ • (t - t )j for St 



< 1 

for St > 1 



(B.4) 



(B.5) 



(B.6) 



for relative velocities due to turbulent motion. 

A comparison between analytical solution and simulation re- 
sult for Brownian motion growth is shown in the top panel of 
Figure A. 2. It can be seen that the position of the peak of the 
grain size distribution (solid line in Figure A. 2) follows the ana- 
lytical result of Eq. B.5. 

A similar comparison for the case of relative velocities due to 
turbulent motion is not as straight-forward since both the turbu- 
lent relative velocities and the dust scale height (Eq. 51) are sub- 
divided into several regimes. We therefore integrated Eq. B.2 nu- 
merically for the case of Brownian motion and turbulent relative 
velocities; the results are shown in the bottom panel of Figure 
A. 2. As before, we see that - after the initial condition is over- 
come - the simulation result follows closely the mono-disperse 
approximation. For grains larger than St = 1, the analytical so- 
lution is also plotted in the bottom panel of Fig. A. 2, almost 
coinciding with the simulation results. 

The radial motion of dust particle was tested in a similar 
fashion: starting from a grain distribution at a radius of 5.5 AU, 
we let the particles drift (taking the gas drag and the radial drift 
into account, see Eq. 19) without coagulation. We compare the 
results to results of a numerical integration of the equation of 
motion for a single particle. The results are shown in Fig. B.l. 

We find that the size distribution behaves as expected: small 
particles are well coupled to the gas, they almost retain their ini- 
tial position since the radial motion due to gas drag is in the 
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order of 0.01 AU. Larger particles, having a larger Stokes num- 
ber drift towards the star on shorter timescales. Particles of a few 



centimeters (corresponding to St=l) are already lost after about 
700 years. 

The mass in all test cases was found to be conserved on the 
order of 10~ n % of the initial value. 



